﻿using System;
using System.Collections.Generic;
using SpiceNet.Diagnostics;
using SpiceNet.Circuits;
using SpiceNet.Simulations;
using SpiceNet.Integration;
using SpiceNet.Parameters;

namespace SpiceNet.Models
{
    /// <summary>
    /// Represents an instance of a bipolar transistor
    /// </summary>
    public class Bipolar : Instance
    {
        /// <summary>
        /// Get parameters
        /// </summary>
        public override Dictionary<string, Parameterized.IP> ParamTable
        {
            get { return _pTable; }
        }

        /// <summary>
        /// Variables that are calculated for the instance
        /// </summary>
        private double BJTtSatCur, BJTtBetaF, BJTtBetaR, BJTtBEleakCur, BJTtBCleakCur, BJTtBEcap, BJTtBEpot, BJTtBCcap, BJTtBCpot, BJTtDepCap, BJTtf1, BJTtf4, BJTtf5, BJTtVcrit;
        private double BJTcapbe, BJTcapbc, BJTcapcs, BJTcapbx;
        private int BJTcolNode, BJTcolPrimeNode, BJTbaseNode, BJTbasePrimeNode, BJTemitNode, BJTemitPrimeNode, BJTsubstNode, BJTstate;
        private Noise.InstanceNoiseData[] BJTnVar = new Noise.InstanceNoiseData[7];

        #region Parameters
        private enum _c
        {
            /* device parameters */
            BJT_AREA = 1,
            BJT_OFF = 2,
            BJT_IC_VBE = 3,
            BJT_IC_VCE = 4,
            BJT_IC = 5,
            BJT_AREA_SENS = 6,
            BJT_TEMP = 7,

            /* device questions */
            BJT_QUEST_FT = 201,
            BJT_QUEST_COLNODE = 202,
            BJT_QUEST_BASENODE = 203,
            BJT_QUEST_EMITNODE = 204,
            BJT_QUEST_SUBSTNODE = 205,
            BJT_QUEST_COLPRIMENODE = 206,
            BJT_QUEST_BASEPRIMENODE = 207,
            BJT_QUEST_EMITPRIMENODE = 208,
            BJT_QUEST_VBE = 209,
            BJT_QUEST_VBC = 210,
            BJT_QUEST_CC = 211,
            BJT_QUEST_CB = 212,
            BJT_QUEST_GPI = 213,
            BJT_QUEST_GMU = 214,
            BJT_QUEST_GM = 215,
            BJT_QUEST_GO = 216,
            BJT_QUEST_QBE = 217,
            BJT_QUEST_CQBE = 218,
            BJT_QUEST_QBC = 219,
            BJT_QUEST_CQBC = 220,
            BJT_QUEST_QCS = 221,
            BJT_QUEST_CQCS = 222,
            BJT_QUEST_QBX = 223,
            BJT_QUEST_CQBX = 224,
            BJT_QUEST_GX = 225,
            BJT_QUEST_CEXBC = 226,
            BJT_QUEST_GEQCB = 227,
            BJT_QUEST_GCCS = 228,
            BJT_QUEST_GEQBX = 229,
            BJT_QUEST_SENS_REAL = 230,
            BJT_QUEST_SENS_IMAG = 231,
            BJT_QUEST_SENS_MAG = 232,
            BJT_QUEST_SENS_PH = 233,
            BJT_QUEST_SENS_CPLX = 234,
            BJT_QUEST_SENS_DC = 235,
            BJT_QUEST_CE = 236,
            BJT_QUEST_CS = 237,
            BJT_QUEST_POWER = 238,
            BJT_QUEST_CPI = 239,
            BJT_QUEST_CMU = 240,
            BJT_QUEST_CBX = 241,
            BJT_QUEST_CCS = 242
        }
        private static Dictionary<string, IP> _pTable = new Dictionary<string, IP>
        {
            { "off", new IP(IF.IOPU, (int)_c.BJT_OFF, T.FLAG, "Device initially off") },
            { "icvbe", new IP(IF.IOPAU, (int)_c.BJT_IC_VBE, T.REAL, "Initial B-E voltage") },
            { "icvce", new IP(IF.IOPAU, (int)_c.BJT_IC_VCE, T.REAL, "Initial C-E voltage") },
            { "area", new IP(IF.IOPU, (int)_c.BJT_AREA, T.REAL, "Area factor") },
            { "init", new IP(IF.IP, (int)_c.BJT_IC, T.REALVEC, "Initial condition vector") }, // We can't use duplicate identifiers... Changed ic to init
            { "sens_area", new IP(IF.IP, (int)_c.BJT_AREA_SENS, T.FLAG, "flag to request sensitivity WRT area") },
            { "colnode", new IP(IF.OPU, (int)_c.BJT_QUEST_COLNODE, T.INTEGER, "Number of collector node") },
            { "basenode", new IP(IF.OPU, (int)_c.BJT_QUEST_BASENODE, T.INTEGER, "Number of base node") },
            { "emitnode", new IP(IF.OPU, (int)_c.BJT_QUEST_EMITNODE, T.INTEGER, "Number of emitter node") },
            { "substnode", new IP(IF.OPU, (int)_c.BJT_QUEST_SUBSTNODE, T.INTEGER, "Number of substrate node") },
            { "colprimenode", new IP(IF.OPU, (int)_c.BJT_QUEST_COLPRIMENODE, T.INTEGER, "Internal collector node") },
            { "baseprimenode", new IP(IF.OPU, (int)_c.BJT_QUEST_BASEPRIMENODE, T.INTEGER, "Internal base node") },
            { "emitprimenode", new IP(IF.OPU, (int)_c.BJT_QUEST_EMITPRIMENODE, T.INTEGER, "Internal emitter node") },
            { "ic", new IP(IF.OP, (int)_c.BJT_QUEST_CC, T.REAL, "Current at collector node") },
            { "ib", new IP(IF.OP, (int)_c.BJT_QUEST_CB, T.REAL, "Current at base node") },
            { "ie", new IP(IF.OP, (int)_c.BJT_QUEST_CE, T.REAL, "Emitter current") },
            { "is", new IP(IF.OPU, (int)_c.BJT_QUEST_CS, T.REAL, "Substrate current") },
            { "vbe", new IP(IF.OP, (int)_c.BJT_QUEST_VBE, T.REAL, "B-E voltage") },
            { "vbc", new IP(IF.OP, (int)_c.BJT_QUEST_VBC, T.REAL, "B-C voltage") },
            { "gm", new IP(IF.OP, (int)_c.BJT_QUEST_GM, T.REAL, "Small signal transconductance") },
            { "gpi", new IP(IF.OP, (int)_c.BJT_QUEST_GPI, T.REAL, "Small signal input conductance - pi") },
            { "gmu", new IP(IF.OP, (int)_c.BJT_QUEST_GMU, T.REAL, "Small signal conductance - mu") },
            { "gx", new IP(IF.OP, (int)_c.BJT_QUEST_GX, T.REAL, "Conductance from base to internal base") },
            { "go", new IP(IF.OP, (int)_c.BJT_QUEST_GO, T.REAL, "Small signal output conductance") },
            { "geqcb", new IP(IF.OPU, (int)_c.BJT_QUEST_GEQCB, T.REAL, "d(Ibe)/d(Vbc)") },
            { "gccs", new IP(IF.OPU, (int)_c.BJT_QUEST_GCCS, T.REAL, "Internal C-S cap. equiv. cond.") },
            { "geqbx", new IP(IF.OPU, (int)_c.BJT_QUEST_GEQBX, T.REAL, "Internal C-B-base cap. equiv. cond.") },
            { "cpi", new IP(IF.OP, (int)_c.BJT_QUEST_CPI, T.REAL, "Internal base to emitter capactance") },
            { "cmu", new IP(IF.OP, (int)_c.BJT_QUEST_CMU, T.REAL, "Internal base to collector capactiance") },
            { "cbx", new IP(IF.OP, (int)_c.BJT_QUEST_CBX, T.REAL, "Base to collector capacitance") },
            { "ccs", new IP(IF.OP, (int)_c.BJT_QUEST_CCS, T.REAL, "Collector to substrate capacitance") },
            { "cqbe", new IP(IF.OPU, (int)_c.BJT_QUEST_CQBE, T.REAL, "Cap. due to charge storage in B-E jct.") },
            { "cqbc", new IP(IF.OPU, (int)_c.BJT_QUEST_CQBC, T.REAL, "Cap. due to charge storage in B-C jct.") },
            { "cqcs", new IP(IF.OPU, (int)_c.BJT_QUEST_CQCS, T.REAL, "Cap. due to charge storage in C-S jct.") },
            { "cqbx", new IP(IF.OPU, (int)_c.BJT_QUEST_CQBX, T.REAL, "Cap. due to charge storage in B-X jct.") },
            { "cexbc", new IP(IF.OPU, (int)_c.BJT_QUEST_CEXBC, T.REAL, "Total Capacitance in B-X junction") },
            { "qbe", new IP(IF.OPU, (int)_c.BJT_QUEST_QBE, T.REAL, "Charge storage B-E junction") },
            { "qbc", new IP(IF.OPU, (int)_c.BJT_QUEST_QBC, T.REAL, "Charge storage B-C junction") },
            { "qcs", new IP(IF.OPU, (int)_c.BJT_QUEST_QCS, T.REAL, "Charge storage C-S junction") },
            { "qbx", new IP(IF.OPU, (int)_c.BJT_QUEST_QBX, T.REAL, "Charge storage B-X junction") },
            { "p", new IP(IF.OPU, (int)_c.BJT_QUEST_POWER, T.REAL, "Power dissipation") },
            { "sens_dc", new IP(IF.OPU, (int)_c.BJT_QUEST_SENS_DC, T.REAL, "dc sensitivity ") },
            { "sens_real", new IP(IF.OPU, (int)_c.BJT_QUEST_SENS_REAL, T.REAL, "real part of ac sensitivity") },
            { "sens_imag", new IP(IF.OPU, (int)_c.BJT_QUEST_SENS_IMAG, T.REAL, "dc sens. & imag part of ac sens.") },
            { "sens_mag", new IP(IF.OPU, (int)_c.BJT_QUEST_SENS_MAG, T.REAL, "sensitivity of ac magnitude") },
            { "sens_ph", new IP(IF.OPU, (int)_c.BJT_QUEST_SENS_PH, T.REAL, "sensitivity of ac phase") },
            { "sens_cplx", new IP(IF.OPU, (int)_c.BJT_QUEST_SENS_CPLX, T.COMPLEX, "ac sensitivity") },
            { "temp", new IP(IF.IOPU, (int)_c.BJT_TEMP, T.REAL, "instance temperature") }
        };

        public Parameter<double> BJTarea { get; } = new Parameter<double>(1.0);
        public Parameter<double> BJTtemp { get; } = new Parameter<double>();
        public bool BJToff { get; set; } = false;
        public Parameter<double> BJTicVBE { get; } = new Parameter<double>();
        public Parameter<double> BJTicVCE { get; } = new Parameter<double>();
        public int BJTsenParmNo { get; set; } = 0;
        #endregion

        #region Constants
        // States
        private const int BJTvbe = 0;
        private const int BJTvbc = 1;
        private const int BJTcc = 2;
        private const int BJTcb = 3;
        private const int BJTgpi = 4;
        private const int BJTgmu = 5;
        private const int BJTgm = 6;
        private const int BJTgo = 7;
        private const int BJTqbe = 8;
        private const int BJTcqbe = 9;
        private const int BJTqbc = 10;
        private const int BJTcqbc = 11;
        private const int BJTqcs = 12;
        private const int BJTcqcs = 13;
        private const int BJTqbx = 14;
        private const int BJTcqbx = 15;
        private const int BJTgx = 16;
        private const int BJTcexbc = 17;
        private const int BJTgeqcb = 18;
        private const int BJTgccs = 19;
        private const int BJTgeqbx = 20;

        private const int BJTRCNOIZ = 0;
        private const int BJTRBNOIZ = 1;
        private const int BJT_RE_NOISE = 2;
        private const int BJTICNOIZ = 3;
        private const int BJTIBNOIZ = 4;
        private const int BJTFLNOIZ = 5;
        private const int BJTTOTNOIZ = 6;
        #endregion

        /// <summary>
        /// Constructor
        /// </summary>
        /// <param name="name"></param>
        public Bipolar(string name) : base(name, 4, typeof(BipolarModel)) { }

        /// <summary>
        /// Constructor
        /// </summary>
        /// <param name="name"></param>
        /// <param name="colnode"></param>
        /// <param name="basenode"></param>
        /// <param name="emitnode"></param>
        /// <param name="substnode"></param>
        public Bipolar(string name, string colnode, string basenode, string emitnode, string substnode) : this(name)
        {
            Connect(colnode, basenode, emitnode, substnode);
        }

        /// <summary>
        /// Set parameter
        /// </summary>
        protected override void Param(int id, object value, Circuit ckt = null)
        {
            double[] v = null;
            switch ((_c)id)
            {
                case _c.BJT_AREA:
                    BJTarea.Par((double)value);
                    break;
                case _c.BJT_TEMP:
                    BJTtemp.Par((double)value + Circuit.CONSTCtoK);
                    break;
                case _c.BJT_OFF:
                    BJToff = (bool)value;
                    break;
                case _c.BJT_IC_VBE:
                    BJTicVBE.Par((double)value);
                    break;
                case _c.BJT_IC_VCE:
                    BJTicVCE.Par((double)value);
                    break;
                case _c.BJT_AREA_SENS:
                    BJTsenParmNo = (int)value;
                    break;
                case _c.BJT_IC:
                    v = (double[])value;
                    switch (v.Length)
                    {
                        case 2:
                            BJTicVCE.Par(v[1]); goto case 1;
                        case 1:
                            BJTicVBE.Par(v[0]); break;
                        default:
                            throw new BadParameterException();
                    }
                    break;
                default:
                    throw new BadParameterException();
            }
        }

        /// <summary>
        /// Ask parameter
        /// </summary>
        protected override object Ask(int id, Circuit ckt = null)
        {
            switch ((_c)id)
            {
                case _c.BJT_QUEST_FT:
                    double tmp = Math.Max(ckt.State.States[0][BJTstate + BJTcqbc],
                        ckt.State.States[0][BJTstate + BJTcqbx]);
                    return BJTgm / (2 * Circuit.CONSTPI *
                        Math.Max(ckt.State.States[0][BJTstate + BJTcqbe], tmp));
                case _c.BJT_TEMP:
                    return BJTtemp - Circuit.CONSTCtoK;
                case _c.BJT_AREA:
                    return BJTarea.Value;
                case _c.BJT_OFF:
                    return BJToff;
                case _c.BJT_IC_VBE:
                    return BJTicVBE.Value;
                case _c.BJT_IC_VCE:
                    return BJTicVCE.Value;
                case _c.BJT_QUEST_COLNODE:
                    return BJTcolNode;
                case _c.BJT_QUEST_BASENODE:
                    return BJTbaseNode;
                case _c.BJT_QUEST_EMITNODE:
                    return BJTemitNode;
                case _c.BJT_QUEST_SUBSTNODE:
                    return BJTsubstNode;
                case _c.BJT_QUEST_COLPRIMENODE:
                    return BJTcolPrimeNode;
                case _c.BJT_QUEST_BASEPRIMENODE:
                    return BJTbasePrimeNode;
                case _c.BJT_QUEST_EMITPRIMENODE:
                    return BJTemitPrimeNode;
                case _c.BJT_QUEST_VBE:
                    return ckt.State.States[0][BJTstate + BJTvbe];
                case _c.BJT_QUEST_VBC:
                    return ckt.State.States[0][BJTstate + BJTvbc];
                case _c.BJT_QUEST_CC:
                    return ckt.State.States[0][BJTstate + BJTcc];
                case _c.BJT_QUEST_CB:
                    return ckt.State.States[0][BJTstate + BJTcb];
                case _c.BJT_QUEST_GPI:
                    return ckt.State.States[0][BJTstate + BJTgpi];
                case _c.BJT_QUEST_GMU:
                    return ckt.State.States[0][BJTstate + BJTgmu];
                case _c.BJT_QUEST_GM:
                    return ckt.State.States[0][BJTstate + BJTgm];
                case _c.BJT_QUEST_GO:
                    return ckt.State.States[0][BJTstate + BJTgo];
                case _c.BJT_QUEST_QBE:
                    return ckt.State.States[0][BJTstate + BJTqbe];
                case _c.BJT_QUEST_CQBE:
                    return ckt.State.States[0][BJTstate + BJTcqbe];
                case _c.BJT_QUEST_QBC:
                    return ckt.State.States[0][BJTstate + BJTqbc];
                case _c.BJT_QUEST_CQBC:
                    return ckt.State.States[0][BJTstate + BJTcqbc];
                case _c.BJT_QUEST_QCS:
                    return ckt.State.States[0][BJTstate + BJTqcs];
                case _c.BJT_QUEST_CQCS:
                    return ckt.State.States[0][BJTstate + BJTcqcs];
                case _c.BJT_QUEST_QBX:
                    return ckt.State.States[0][BJTstate + BJTqbx];
                case _c.BJT_QUEST_CQBX:
                    return ckt.State.States[0][BJTstate + BJTcqbx];
                case _c.BJT_QUEST_GX:
                    return ckt.State.States[0][BJTstate + BJTgx];
                case _c.BJT_QUEST_CEXBC:
                    return ckt.State.States[0][BJTstate + BJTcexbc];
                case _c.BJT_QUEST_GEQCB:
                    return ckt.State.States[0][BJTstate + BJTgeqcb];
                case _c.BJT_QUEST_GCCS:
                    return ckt.State.States[0][BJTstate + BJTgccs];
                case _c.BJT_QUEST_GEQBX:
                    return ckt.State.States[0][BJTstate + BJTgeqbx];
                case _c.BJT_QUEST_CS:
                    if (ckt.CurrentAnalysis.HasFlag(AnalysisType.DoingAc))
                    {
                        throw new AcCurrentPowerException();
                    }
                    else if ((ckt.CurrentAnalysis & (AnalysisType.DoingDcOp | AnalysisType.DoingTrCv)) != 0)
                    {
                        return 0;
                    }
                    else if ((ckt.CurrentAnalysis.HasFlag(AnalysisType.DoingTran)) &&
                             (ckt.Mode.HasFlag(Circuit.Modes.TranOp)))
                    {
                        return 0;
                    }
                    else
                    {
                        return -ckt.State.States[0][BJTstate + BJTcqcs];
                    }
                case _c.BJT_QUEST_CE:
                    if (ckt.CurrentAnalysis.HasFlag(AnalysisType.DoingAc))
                    {
                        throw new AcCurrentPowerException();
                    }
                    else
                    {
                        double value = -ckt.State.States[0][BJTstate + BJTcc];
                        value -= ckt.State.States[0][BJTstate + BJTcb];
                        if ((ckt.CurrentAnalysis.HasFlag(AnalysisType.DoingTran)) &&
                                !(ckt.Mode.HasFlag(Circuit.Modes.TranOp)))
                        {
                            value += ckt.State.States[0][BJTstate + BJTcqcs];
                        }
                        return value;
                    }
                case _c.BJT_QUEST_POWER:
                    if (ckt.CurrentAnalysis.HasFlag(AnalysisType.DoingAc))
                    {
                        throw new AcCurrentPowerException();
                    }
                    else
                    {
                        double value = ckt.State.States[0][BJTstate + BJTcc] *
                                ckt.State.Solution[BJTcolNode];
                        value += ckt.State.States[0][BJTstate + BJTcb] *
                                ckt.State.Solution[BJTbaseNode];
                        if ((ckt.CurrentAnalysis.HasFlag(AnalysisType.DoingTran)) && !(ckt.Mode.HasFlag(Circuit.Modes.TranOp)))
                        {
                            value -= ckt.State.States[0][BJTstate + BJTcqcs] *
                                    ckt.State.Solution[BJTsubstNode];
                        }
                        tmp = -ckt.State.States[0][BJTstate + BJTcc];
                        tmp -= ckt.State.States[0][BJTstate + BJTcb];
                        if ((ckt.CurrentAnalysis.HasFlag(AnalysisType.DoingTran)) &&
                                !(ckt.Mode.HasFlag(Circuit.Modes.TranOp)))
                        {
                            tmp += ckt.State.States[0][BJTstate + BJTcqcs];
                        }
                        value += tmp * ckt.State.Solution[BJTemitNode];
                        return value;
                    }
                case _c.BJT_QUEST_CPI:
                    return BJTcapbe;
                case _c.BJT_QUEST_CMU:
                    return BJTcapbc;
                case _c.BJT_QUEST_CBX:
                    return BJTcapbx;
                case _c.BJT_QUEST_CCS:
                    return BJTcapcs;
                default:
                    throw new BadParameterException();
            }
        }

        /// <summary>
        /// Set up the bipolar transistor
        /// </summary>
        /// <param name="ckt"></param>
        public override void Setup(Circuit ckt, Model mod, object p = null)
        {
            BipolarModel model = mod as BipolarModel;

            // Add extra nodes if the bipolar transistor has conductances at the collector or emitter
            List<NodeType> extra = new List<NodeType>();
            if (model.BJTcollectorResist > 0.0)
                extra.Add(NodeType.Voltage);
            if (model.BJTemitterResist > 0.0)
                extra.Add(NodeType.Voltage);
            if (model.BJTbaseResist > 0.0)
                extra.Add(NodeType.Voltage);

            // Bind the nodes
            CircuitNode[] nodes = BindNodes(ckt, extra.ToArray());
            int index = 4;
            BJTcolNode = nodes[0].Number;
            BJTbaseNode = nodes[1].Number;
            BJTemitNode = nodes[2].Number;
            BJTsubstNode = nodes[3].Number;
            BJTcolPrimeNode = model.BJTcollectorResist == 0.0 ? BJTcolNode : nodes[index++].Number;
            BJTbasePrimeNode = model.BJTbaseResist == 0.0 ? BJTbaseNode : nodes[index++].Number;
            BJTemitPrimeNode = model.BJTemitterResist == 0.0 ? BJTemitNode : nodes[index++].Number;

            // Reserve states
            BJTstate = CreateState(ckt, 21);

            // TODO: Sensitivity
        }

        /// <summary>
        /// Do temperature-dependent stuff
        /// </summary>
        /// <param name="ckt"></param>
        public override void Temperature(Circuit ckt, Model mod, object p = null)
        {
            BipolarModel model = mod as BipolarModel;
            var tp = p as BipolarModel.TempParams;
            double vt;
            double ratlog;
            double ratio1;
            double factlog;
            double bfactor;
            double factor;
            double fact2;
            double pbo;
            double pbfact;
            double gmaold;
            double gmanew;
            double egfet;
            double arg;

            if (!BJTtemp.Given) BJTtemp.Value = ckt.Config.Temperature;
            vt = BJTtemp * Circuit.CONSTKoverQ;
            fact2 = BJTtemp / Circuit.CONSTRefTemp;
            egfet = 1.16 - (7.02e-4 * BJTtemp * BJTtemp) /
                    (BJTtemp + 1108);
            arg = -egfet / (2 * Circuit.CONSTBoltz * BJTtemp) +
                    1.1150877 / (Circuit.CONSTBoltz * (Circuit.CONSTRefTemp + Circuit.CONSTRefTemp));
            pbfact = -2 * vt * (1.5 * Math.Log(fact2) + Circuit.CHARGE * arg);

            ratlog = Math.Log(BJTtemp / model.BJTtnom);
            ratio1 = BJTtemp / model.BJTtnom - 1;
            factlog = ratio1 * model.BJTenergyGap / vt +
                    model.BJTtempExpIS * ratlog;
            factor = Math.Exp(factlog);
            BJTtSatCur = model.BJTsatCur * factor;
            bfactor = Math.Exp(ratlog * model.BJTbetaExp);
            BJTtBetaF = model.BJTbetaF * bfactor;
            BJTtBetaR = model.BJTbetaR * bfactor;
            BJTtBEleakCur = model.BJTleakBEcurrent *
                    Math.Exp(factlog / model.BJTleakBEemissionCoeff) / bfactor;
            BJTtBCleakCur = model.BJTleakBCcurrent *
                    Math.Exp(factlog / model.BJTleakBCemissionCoeff) / bfactor;

            pbo = (model.BJTpotentialBE - pbfact) / tp.fact1;
            gmaold = (model.BJTpotentialBE - pbo) / pbo;
            BJTtBEcap = model.BJTdepletionCapBE /
                    (1 + model.BJTjunctionExpBE *
                    (4e-4 * (model.BJTtnom - Circuit.CONSTRefTemp) - gmaold));
            BJTtBEpot = fact2 * pbo + pbfact;
            gmanew = (BJTtBEpot - pbo) / pbo;
            BJTtBEcap *= 1 + model.BJTjunctionExpBE *
                    (4e-4 * (BJTtemp - Circuit.CONSTRefTemp) - gmanew);

            pbo = (model.BJTpotentialBC - pbfact) / tp.fact1;
            gmaold = (model.BJTpotentialBC - pbo) / pbo;
            BJTtBCcap = model.BJTdepletionCapBC /
                    (1 + model.BJTjunctionExpBC *
                    (4e-4 * (model.BJTtnom - Circuit.CONSTRefTemp) - gmaold));
            BJTtBCpot = fact2 * pbo + pbfact;
            gmanew = (BJTtBCpot - pbo) / pbo;
            BJTtBCcap *= 1 + model.BJTjunctionExpBC *
                    (4e-4 * (BJTtemp - Circuit.CONSTRefTemp) - gmanew);

            BJTtDepCap = model.BJTdepletionCapCoeff * BJTtBEpot;
            BJTtf1 = BJTtBEpot * (1 - Math.Exp((1 -
                    model.BJTjunctionExpBE) * tp.xfc)) /
                    (1 - model.BJTjunctionExpBE);
            BJTtf4 = model.BJTdepletionCapCoeff * BJTtBCpot;
            BJTtf5 = BJTtBCpot * (1 - Math.Exp((1 -
                    model.BJTjunctionExpBC) * tp.xfc)) /
                    (1 - model.BJTjunctionExpBC);
            BJTtVcrit = vt * Math.Log(vt / (Circuit.CONSTroot2 * model.BJTsatCur));
        }

        /// <summary>
        /// Load the BJT transistor in the circuit
        /// </summary>
        /// <param name="ckt"></param>
        public override void Load(Circuit ckt, Model mod, object p = null)
        {
            BipolarModel model = mod as BipolarModel;
            CircuitState state = ckt.State;
            var integ = ckt.Integration;
            IntegrationMethod.IntegrationResult result;

            double arg1;
            double arg2;
            double arg3;
            double arg;
            double argtf;
            double c2;
            double c4;
            double capbc;
            double capbe;
            double cb;
            double cbc;
            double cbcn;
            double cbe;
            double cben;
            double cbhat;
            double cc;
            double cchat;
            double cdis;
            double ceq;
            double ceqbc;
            double ceqbe;
            double ceqbx;
            double ceqcs;
            double cex;
            double csat;
            double ctot;
            double czbc;
            double czbcf2;
            double czbe;
            double czbef2;
            double czbx;
            double czbxf2;
            double czcs;
            double delvbc;
            double delvbe;
            double denom;
            double dqbdvc;
            double dqbdve;
            double evbc;
            double evbcn;
            double evbe;
            double evben;
            double f1;
            double f2;
            double f3;
            double fcpc;
            double fcpe;
            double gbc;
            double gbcn;
            double gbe;
            double gben;
            double gccs;
            double gcpr;
            double gepr;
            double geq;
            double geqbx;
            double geqcb;
            double gex;
            double gm;
            double gmu;
            double go;
            double gpi;
            double gx;
            double oik;
            double oikr;
            double ovtf;
            double pc;
            double pe;
            double ps;
            double q1;
            double q2;
            double qb;
            double rbpi;
            double rbpr;
            double sarg;
            double sqarg;
            double td;
            double temp;
            double tf;
            double tr;
            double vbc;
            double vbe;
            double vbx;
            double vce;
            double vcs;
            double vt;
            double vtc;
            double vte;
            double vtn;
            double xjrb;
            double xjtf;
            double xmc;
            double xme;
            double xms;
            double xtf;
            double capbx = 0;
            double capcs = 0;
            bool icheck, ichk1;


            vt = BJTtemp * Circuit.CONSTKoverQ;

            gccs = 0;
            ceqcs = 0;
            geqbx = 0;
            ceqbx = 0;
            geqcb = 0;
            /*
             *   dc model paramters
             */
            csat = BJTtSatCur * BJTarea;
            rbpr = model.BJTminBaseResist / BJTarea;
            rbpi = model.BJTbaseResist / BJTarea - rbpr;
            gcpr = model.BJTcollectorConduct * BJTarea;
            gepr = model.BJTemitterConduct * BJTarea;
            oik = model.BJTinvRollOffF / BJTarea;
            c2 = BJTtBEleakCur * BJTarea;
            vte = model.BJTleakBEemissionCoeff * vt;
            oikr = model.BJTinvRollOffR / BJTarea;
            c4 = BJTtBCleakCur * BJTarea;
            vtc = model.BJTleakBCemissionCoeff * vt;
            td = model.BJTexcessPhaseFactor;
            xjrb = model.BJTbaseCurrentHalfResist * BJTarea;

            /*
             *   initialization
             */
            icheck = true;
            if (ckt.Mode.HasFlag(Circuit.Modes.InitSmSig))
            {
                vbe = state.States[0][BJTstate + BJTvbe];
                vbc = state.States[0][BJTstate + BJTvbc];
                vbx = model.BJTtype * (
                    state.OldSolution[BJTbaseNode] -
                    state.OldSolution[BJTcolPrimeNode]);
                vcs = model.BJTtype * (
                    state.OldSolution[BJTsubstNode] -
                    state.OldSolution[BJTcolPrimeNode]);
            }
            else if (ckt.Mode.HasFlag(Circuit.Modes.InitTran))
            {
                vbe = state.States[1][BJTstate + BJTvbe];
                vbc = state.States[1][BJTstate + BJTvbc];
                vbx = model.BJTtype * (
                    state.OldSolution[BJTbaseNode] -
                    state.OldSolution[BJTcolPrimeNode]);
                vcs = model.BJTtype * (
                    state.OldSolution[BJTsubstNode] -
                    state.OldSolution[BJTcolPrimeNode]);
                if ((ckt.Mode.HasFlag(Circuit.Modes.Tran)) && (ckt.Mode.HasFlag(Circuit.Modes.Uic)))
                {
                    vbx = model.BJTtype * (BJTicVBE - BJTicVCE);
                    vcs = 0;
                }
            }
            else if ((ckt.Mode.HasFlag(Circuit.Modes.InitJct)) &&
                  (ckt.Mode.HasFlag(Circuit.Modes.TranOp)) && (ckt.Mode.HasFlag(Circuit.Modes.Uic)))
            {
                vbe = model.BJTtype * BJTicVBE;
                vce = model.BJTtype * BJTicVCE;
                vbc = vbe - vce;
                vbx = vbc;
                vcs = 0;
            }
            else if ((ckt.Mode.HasFlag(Circuit.Modes.InitJct)) && (!BJToff))
            {
                vbe = BJTtVcrit;
                vbc = 0;
                /* ERROR:  need to initialize VCS, VBX here */
                vcs = vbx = 0;
            }
            else if ((ckt.Mode.HasFlag(Circuit.Modes.InitJct)) ||
                  ((ckt.Mode.HasFlag(Circuit.Modes.InitFix)) && (BJToff)))
            {
                vbe = 0;
                vbc = 0;
                /* ERROR:  need to initialize VCS, VBX here */
                vcs = vbx = 0;
            }
            else
            {
                /*
                 *   compute new nonlinear branch voltages
                 */
                vbe = model.BJTtype * (
                    state.OldSolution[BJTbasePrimeNode] -
                    state.OldSolution[BJTemitPrimeNode]);
                vbc = model.BJTtype * (
                    state.OldSolution[BJTbasePrimeNode] -
                    state.OldSolution[BJTcolPrimeNode]);
                delvbe = vbe - state.States[0][BJTstate + BJTvbe];
                delvbc = vbc - state.States[0][BJTstate + BJTvbc];
                vbx = model.BJTtype * (
                    state.OldSolution[BJTbaseNode] -
                    state.OldSolution[BJTcolPrimeNode]);
                vcs = model.BJTtype * (
                    state.OldSolution[BJTsubstNode] -
                    state.OldSolution[BJTcolPrimeNode]);
                cchat = state.States[0][BJTstate + BJTcc] + (state.States[0][BJTstate + BJTgm] + state.States[0][BJTstate + BJTgo]) * delvbe -
                        (state.States[0][BJTstate + BJTgo] + state.States[0][BJTstate + BJTgmu]) * delvbc;
                cbhat = state.States[0][BJTstate + BJTcb] + state.States[0][BJTstate + BJTgpi] * delvbe + state.States[0][BJTstate + BJTgmu] *
                        delvbc;
                /*
                 *   limit nonlinear branch voltages
                 */
                ichk1 = true;
                vbe = pnjlim(vbe, state.States[0][BJTstate + BJTvbe], vt,
                        BJTtVcrit, ref icheck);
                vbc = pnjlim(vbc, state.States[0][BJTstate + BJTvbc], vt,
                        BJTtVcrit, ref ichk1);
                if (ichk1) icheck = true;
            }
        /*
         *   determine dc current and derivitives
         */
        // next1: 
            vtn = vt * model.BJTemissionCoeffF;
            if (vbe > -5 * vtn)
            {
                evbe = Math.Exp(vbe / vtn);
                cbe = csat * (evbe - 1) + ckt.Config.Gmin * vbe;
                gbe = csat * evbe / vtn + ckt.Config.Gmin;
                if (c2 == 0)
                {
                    cben = 0;
                    gben = 0;
                }
                else
                {
                    evben = Math.Exp(vbe / vte);
                    cben = c2 * (evben - 1);
                    gben = c2 * evben / vte;
                }
            }
            else
            {
                gbe = -csat / vbe + ckt.Config.Gmin;
                cbe = gbe * vbe;
                gben = -c2 / vbe;
                cben = gben * vbe;
            }
            vtn = vt * model.BJTemissionCoeffR;
            if (vbc > -5 * vtn)
            {
                evbc = Math.Exp(vbc / vtn);
                cbc = csat * (evbc - 1) + ckt.Config.Gmin * vbc;
                gbc = csat * evbc / vtn + ckt.Config.Gmin;
                if (c4 == 0)
                {
                    cbcn = 0;
                    gbcn = 0;
                }
                else
                {
                    evbcn = Math.Exp(vbc / vtc);
                    cbcn = c4 * (evbcn - 1);
                    gbcn = c4 * evbcn / vtc;
                }
            }
            else
            {
                gbc = -csat / vbc + ckt.Config.Gmin;
                cbc = gbc * vbc;
                gbcn = -c4 / vbc;
                cbcn = gbcn * vbc;
            }
            /*
             *   determine base charge terms
             */
            q1 = 1 / (1 - model.BJTinvEarlyVoltF * vbc - model.BJTinvEarlyVoltR * vbe);
            if (oik == 0 && oikr == 0)
            {
                qb = q1;
                dqbdve = q1 * qb * model.BJTinvEarlyVoltR;
                dqbdvc = q1 * qb * model.BJTinvEarlyVoltF;
            }
            else
            {
                q2 = oik * cbe + oikr * cbc;
                arg = Math.Max(0, 1 + 4 * q2);
                sqarg = 1;
                if (arg != 0) sqarg = Math.Sqrt(arg);
                qb = q1 * (1 + sqarg) / 2;
                dqbdve = q1 * (qb * model.BJTinvEarlyVoltR + oik * gbe / sqarg);
                dqbdvc = q1 * (qb * model.BJTinvEarlyVoltF + oikr * gbc / sqarg);
            }
            /*
             *   weil's approx. for excess phase applied with backward-
             *   euler integration
             */
            cc = 0;
            cex = cbe;
            gex = gbe;
            if ((ckt.Mode & (Circuit.Modes.Tran | Circuit.Modes.Ac)) != 0 && td != 0)
            {
                arg1 = integ.Delta / td;
                arg2 = 3 * arg1;
                arg1 = arg2 * arg1;
                denom = 1 + arg1 + arg2;
                arg3 = arg1 / denom;
                if (ckt.Mode.HasFlag(Circuit.Modes.InitTran))
                {
                    state.States[1][BJTstate + BJTcexbc] = cbe / qb;
                    state.States[2][BJTstate + BJTcexbc] =
                            state.States[1][BJTstate + BJTcexbc];
                }
                cc = (state.States[1][BJTstate + BJTcexbc] * (1 + integ.Delta /
                        integ.DeltaOld[1] + arg2) -
                        state.States[2][BJTstate + BJTcexbc] * integ.Delta /
                        integ.DeltaOld[1]) / denom;
                cex = cbe * arg3;
                gex = gbe * arg3;
                state.States[0][BJTstate + BJTcexbc] = cc + cex / qb;
            }
            /*
             *   determine dc incremental conductances
             */
            cc = cc + (cex - cbc) / qb - cbc / BJTtBetaR - cbcn;
            cb = cbe / BJTtBetaF + cben + cbc / BJTtBetaR + cbcn;
            gx = rbpr + rbpi / qb;
            if (xjrb != 0)
            {
                arg1 = Math.Max(cb / xjrb, 1e-9);
                arg2 = (-1 + Math.Sqrt(1 + 14.59025 * arg1)) / 2.4317 / Math.Sqrt(arg1);
                arg1 = Math.Tan(arg2);
                gx = rbpr + 3 * rbpi * (arg1 - arg2) / arg2 / arg1 / arg1;
            }
            if (gx != 0) gx = 1 / gx;
            gpi = gbe / BJTtBetaF + gben;
            gmu = gbc / BJTtBetaR + gbcn;
            go = (gbc + (cex - cbc) * dqbdvc / qb) / qb;
            gm = (gex - (cex - cbc) * dqbdve / qb) / qb - go;
            if (((ckt.Mode & (Circuit.Modes.Tran | Circuit.Modes.Ac)) != 0) ||
                    ((ckt.Mode.HasFlag(Circuit.Modes.TranOp)) && (ckt.Mode.HasFlag(Circuit.Modes.Uic))) ||
                    (ckt.Mode.HasFlag(Circuit.Modes.InitSmSig)))
            {
                /*
                 *   charge storage elements
                 */
                tf = model.BJTtransitTimeF;
                tr = model.BJTtransitTimeR;
                czbe = BJTtBEcap * BJTarea;
                pe = BJTtBEpot;
                xme = model.BJTjunctionExpBE;
                cdis = model.BJTbaseFractionBCcap;
                ctot = BJTtBCcap * BJTarea;
                czbc = ctot * cdis;
                czbx = ctot - czbc;
                pc = BJTtBCpot;
                xmc = model.BJTjunctionExpBC;
                fcpe = BJTtDepCap;
                czcs = model.BJTcapCS * BJTarea;
                ps = model.BJTpotentialSubstrate;
                xms = model.BJTexponentialSubstrate;
                xtf = model.BJTtransitTimeBiasCoeffF;
                ovtf = model.BJTtransitTimeVBCFactor;
                xjtf = model.BJTtransitTimeHighCurrentF * BJTarea;
                if (tf != 0 && vbe > 0)
                {
                    argtf = 0;
                    arg2 = 0;
                    arg3 = 0;
                    if (xtf != 0)
                    {
                        argtf = xtf;
                        if (ovtf != 0)
                        {
                            argtf = argtf * Math.Exp(vbc * ovtf);
                        }
                        arg2 = argtf;
                        if (xjtf != 0)
                        {
                            temp = cbe / (cbe + xjtf);
                            argtf = argtf * temp * temp;
                            arg2 = argtf * (3 - temp - temp);
                        }
                        arg3 = cbe * argtf * ovtf;
                    }
                    cbe = cbe * (1 + argtf) / qb;
                    gbe = (gbe * (1 + arg2) - cbe * dqbdve) / qb;
                    geqcb = tf * (arg3 - cbe * dqbdvc) / qb;
                }
                if (vbe < fcpe)
                {
                    arg = 1 - vbe / pe;
                    sarg = Math.Exp(-xme * Math.Log(arg));
                    state.States[0][BJTstate + BJTqbe] = tf * cbe + pe * czbe *
                            (1 - arg * sarg) / (1 - xme);
                    capbe = tf * gbe + czbe * sarg;
                }
                else
                {
                    f1 = BJTtf1;
                    f2 = model.BJTf2;
                    f3 = model.BJTf3;
                    czbef2 = czbe / f2;
                    state.States[0][BJTstate + BJTqbe] = tf * cbe + czbe * f1 + czbef2 *
                            (f3 * (vbe - fcpe) + (xme / (pe + pe)) * (vbe * vbe - fcpe * fcpe));
                    capbe = tf * gbe + czbef2 * (f3 + xme * vbe / pe);
                }
                fcpc = BJTtf4;
                f1 = BJTtf5;
                f2 = model.BJTf6;
                f3 = model.BJTf7;
                if (vbc < fcpc)
                {
                    arg = 1 - vbc / pc;
                    sarg = Math.Exp(-xmc * Math.Log(arg));
                    state.States[0][BJTstate + BJTqbc] = tr * cbc + pc * czbc * (
                            1 - arg * sarg) / (1 - xmc);
                    capbc = tr * gbc + czbc * sarg;
                }
                else
                {
                    czbcf2 = czbc / f2;
                    state.States[0][BJTstate + BJTqbc] = tr * cbc + czbc * f1 + czbcf2 *
                            (f3 * (vbc - fcpc) + (xmc / (pc + pc)) * (vbc * vbc - fcpc * fcpc));
                    capbc = tr * gbc + czbcf2 * (f3 + xmc * vbc / pc);
                }
                if (vbx < fcpc)
                {
                    arg = 1 - vbx / pc;
                    sarg = Math.Exp(-xmc * Math.Log(arg));
                    state.States[0][BJTstate + BJTqbx] =
                        pc * czbx * (1 - arg * sarg) / (1 - xmc);
                    capbx = czbx * sarg;
                }
                else
                {
                    czbxf2 = czbx / f2;
                    state.States[0][BJTstate + BJTqbx] = czbx * f1 + czbxf2 *
                            (f3 * (vbx - fcpc) + (xmc / (pc + pc)) * (vbx * vbx - fcpc * fcpc));
                    capbx = czbxf2 * (f3 + xmc * vbx / pc);
                }
                if (vcs < 0)
                {
                    arg = 1 - vcs / ps;
                    sarg = Math.Exp(-xms * Math.Log(arg));
                    state.States[0][BJTstate + BJTqcs] = ps * czcs * (1 - arg * sarg) /
                            (1 - xms);
                    capcs = czcs * sarg;
                }
                else
                {
                    state.States[0][BJTstate + BJTqcs] = vcs * czcs * (1 + xms * vcs /
                            (2 * ps));
                    capcs = czcs * (1 + xms * vcs / ps);
                }
                BJTcapbe = capbe;
                BJTcapbc = capbc;
                BJTcapcs = capcs;
                BJTcapbx = capbx;

                /*
                 *   store small-signal parameters
                 */
                if ((!(ckt.Mode.HasFlag(Circuit.Modes.TranOp))) ||
                        (!(ckt.Mode.HasFlag(Circuit.Modes.Uic))))
                {
                    if (ckt.Mode.HasFlag(Circuit.Modes.InitSmSig))
                    {
                        state.States[0][BJTstate + BJTcqbe] = capbe;
                        state.States[0][BJTstate + BJTcqbc] = capbc;
                        state.States[0][BJTstate + BJTcqcs] = capcs;
                        state.States[0][BJTstate + BJTcqbx] = capbx;
                        state.States[0][BJTstate + BJTcexbc] = geqcb;
                        return; /* go to 1000 */
                    }
                    /*
                     *   transient analysis
                     */

                    if (ckt.Mode.HasFlag(Circuit.Modes.InitTran))
                    {
                        state.States[1][BJTstate + BJTqbe] =
                                state.States[0][BJTstate + BJTqbe];
                        state.States[1][BJTstate + BJTqbc] =
                                state.States[0][BJTstate + BJTqbc];
                        state.States[1][BJTstate + BJTqbx] =
                                state.States[0][BJTstate + BJTqbx];
                        state.States[1][BJTstate + BJTqcs] =
                                state.States[0][BJTstate + BJTqcs];
                    }
                    result = integ.Integrate(state, BJTstate + BJTqbe, capbe);
                    geq = result.Geq;
                    ceq = result.Ceq;
                    geqcb = geqcb * integ.Deriv;
                    gpi = gpi + geq;
                    cb = cb + state.States[0][BJTstate + BJTcqbe];
                    result = integ.Integrate(state, BJTstate + BJTqbc, capbc);
                    geq = result.Geq;
                    ceq = result.Ceq;
                    gmu = gmu + geq;
                    cb = cb + state.States[0][BJTstate + BJTcqbc];
                    cc = cc - state.States[0][BJTstate + BJTcqbc];
                    if (ckt.Mode.HasFlag(Circuit.Modes.InitTran))
                    {
                        state.States[1][BJTstate + BJTcqbe] =
                                state.States[0][BJTstate + BJTcqbe];
                        state.States[1][BJTstate + BJTcqbc] =
                                state.States[0][BJTstate + BJTcqbc];
                    }
                }
            }

            /*
             *   check convergence
             */
            if ((!(ckt.Mode.HasFlag(Circuit.Modes.InitFix))) || (!(BJToff)))
            {
                if (icheck)
                {
                    ckt.IsCon = false;
                }
            }

            /*
             *      charge storage for c-s and b-x junctions
             */
            if ((ckt.Mode & (Circuit.Modes.Tran | Circuit.Modes.Ac)) != 0)
            {
                result = integ.Integrate(state, BJTstate + BJTqcs, capcs);
                gccs = result.Geq;
                ceq = result.Ceq;
                result = integ.Integrate(state, BJTstate + BJTqbx, capbx);
                geqbx = result.Geq;
                ceq = result.Ceq;
                if (ckt.Mode.HasFlag(Circuit.Modes.InitTran))
                {
                    state.States[1][BJTstate + BJTcqbx] =
                            state.States[0][BJTstate + BJTcqbx];
                    state.States[1][BJTstate + BJTcqcs] =
                            state.States[0][BJTstate + BJTcqcs];
                }
            }
        // next2:
            state.States[0][BJTstate + BJTvbe] = vbe;
            state.States[0][BJTstate + BJTvbc] = vbc;
            state.States[0][BJTstate + BJTcc] = cc;
            state.States[0][BJTstate + BJTcb] = cb;
            state.States[0][BJTstate + BJTgpi] = gpi;
            state.States[0][BJTstate + BJTgmu] = gmu;
            state.States[0][BJTstate + BJTgm] = gm;
            state.States[0][BJTstate + BJTgo] = go;
            state.States[0][BJTstate + BJTgx] = gx;
            state.States[0][BJTstate + BJTgeqcb] = geqcb;
            state.States[0][BJTstate + BJTgccs] = gccs;
            state.States[0][BJTstate + BJTgeqbx] = geqbx;

            /* Do not load the Jacobian and the rhs if
               perturbation is being carried out */

        // load:
            /*
             *  load current excitation vector
             */
            ceqcs = model.BJTtype * (state.States[0][BJTstate + BJTcqcs] -
                    vcs * gccs);
            ceqbx = model.BJTtype * (state.States[0][BJTstate + BJTcqbx] -
                    vbx * geqbx);
            ceqbe = model.BJTtype * (cc + cb - vbe * (gm + go + gpi) + vbc *
                    (go - geqcb));
            ceqbc = model.BJTtype * (-cc + vbe * (gm + go) - vbc * (gmu + go));

            state.Rhs[BJTbaseNode] += (-ceqbx);
            state.Rhs[BJTcolPrimeNode] +=
                    (ceqcs + ceqbx + ceqbc);
            state.Rhs[BJTbasePrimeNode] +=
                    (-ceqbe - ceqbc);
            state.Rhs[BJTemitPrimeNode] += (ceqbe);
            state.Rhs[BJTsubstNode] += (-ceqcs);
            /*
             *  load y matrix
             */
            state.Matrix[BJTcolNode, BJTcolNode] += (gcpr);
            state.Matrix[BJTbaseNode, BJTbaseNode] += (gx + geqbx);
            state.Matrix[BJTemitNode, BJTemitNode] += (gepr);
            state.Matrix[BJTcolPrimeNode, BJTcolPrimeNode] += (gmu + go + gcpr + gccs + geqbx);
            state.Matrix[BJTbasePrimeNode, BJTbasePrimeNode] += (gx + gpi + gmu + geqcb);
            state.Matrix[BJTemitPrimeNode, BJTemitPrimeNode] += (gpi + gepr + gm + go);
            state.Matrix[BJTcolNode, BJTcolPrimeNode] += (-gcpr);
            state.Matrix[BJTbaseNode, BJTbasePrimeNode] += (-gx);
            state.Matrix[BJTemitNode, BJTemitPrimeNode] += (-gepr);
            state.Matrix[BJTcolPrimeNode, BJTcolNode] += (-gcpr);
            state.Matrix[BJTcolPrimeNode, BJTbasePrimeNode] += (-gmu + gm);
            state.Matrix[BJTcolPrimeNode, BJTemitPrimeNode] += (-gm - go);
            state.Matrix[BJTbasePrimeNode, BJTbaseNode] += (-gx);
            state.Matrix[BJTbasePrimeNode, BJTcolPrimeNode] += (-gmu - geqcb);
            state.Matrix[BJTbasePrimeNode, BJTemitPrimeNode] += (-gpi);
            state.Matrix[BJTemitPrimeNode, BJTemitNode] += (-gepr);
            state.Matrix[BJTemitPrimeNode, BJTcolPrimeNode] += (-go + geqcb);
            state.Matrix[BJTemitPrimeNode, BJTbasePrimeNode] += (-gpi - gm - geqcb);
            state.Matrix[BJTsubstNode, BJTsubstNode] += (gccs);
            state.Matrix[BJTcolPrimeNode, BJTsubstNode] += (-gccs);
            state.Matrix[BJTsubstNode, BJTcolPrimeNode] += (-gccs);
            state.Matrix[BJTbaseNode, BJTcolPrimeNode] += (-geqbx);
            state.Matrix[BJTcolPrimeNode, BJTbaseNode] += (-geqbx);
        }

        /// <summary>
        /// Load the BJT transistor in the circuit for AC or Pole-Zero analysis
        /// </summary>
        /// <param name="ckt"></param>
        public override void AcLoad(Circuit ckt, Model mod, object p = null)
        {
            BipolarModel model = mod as BipolarModel;

            CircuitStateComplex state = ckt.StateComplex;
            CircuitState cstate = ckt.State;

            double gcpr;
            double gepr;
            double gpi;
            double gmu;
            double go;
            double xgm;
            // double td;
            // double arg;
            double gm;
            double gx;
            double xcpi;
            double xcmu;
            double xcbx;
            double xccs;
            double xcmcb;

            gcpr = model.BJTcollectorConduct * BJTarea; // [EDIT] Changed model.BJTcollectorResist to model.BJTcollectorConduct
            gepr = model.BJTemitterConduct * BJTarea; // [EDIT] Changed model.BJTemitterResist to model.BJTemitterConduct
            gpi = cstate.States[0][BJTstate + BJTgpi];
            gmu = cstate.States[0][BJTstate + BJTgmu];
            gm = cstate.States[0][BJTstate + BJTgm];
            go = cstate.States[0][BJTstate + BJTgo];
            xgm = 0;
            gx = cstate.States[0][BJTstate + BJTgx];
            xcpi = cstate.States[0][BJTstate + BJTcqbe];
            xcmu = cstate.States[0][BJTstate + BJTcqbc];
            xcbx = cstate.States[0][BJTstate + BJTcqbx];
            xccs = cstate.States[0][BJTstate + BJTcqcs];
            xcmcb = cstate.States[0][BJTstate + BJTcexbc];

            state.Matrix[BJTcolNode, BJTcolNode] += gcpr;
            // *(BJTcolColPtr) += (gcpr);
            state.Matrix[BJTbaseNode, BJTbaseNode] += gx + xcbx * state.Omega;
            // *(BJTbaseBasePtr) += (gx) + (xcbx) * (s->real);
            // *(BJTbaseBasePtr + 1) += (xcbx) * (s->imag);

            state.Matrix[BJTemitNode, BJTemitNode] += gepr;
            // *(BJTemitEmitPtr) += (gepr);
            state.Matrix[BJTcolPrimeNode, BJTcolPrimeNode] += (gmu + go + gcpr) + (xcmu + xccs + xcbx) * state.Omega;
            // *(BJTcolPrimeColPrimePtr) += (gmu + go + gcpr)
            //                             + (xcmu + xccs + xcbx) * (s->real);
            // *(BJTcolPrimeColPrimePtr + 1) += (xcmu + xccs + xcbx)
            //                                      * (s->imag);
            state.Matrix[BJTbasePrimeNode, BJTbasePrimeNode] += (gx + gpi + gmu) + (xcpi + xcmu + xcmcb) * state.Omega;
            // *(BJTbasePrimeBasePrimePtr) += (gx + gpi + gmu)
            //                        + (xcpi + xcmu + xcmcb) * (s->real);
            // *(BJTbasePrimeBasePrimePtr + 1) += (xcpi + xcmu + xcmcb)
            //                                        * (s->imag);
            state.Matrix[BJTemitPrimeNode, BJTemitPrimeNode] += (gpi + gepr + gm + go) + (xcpi + xgm) * state.Omega;
            // *(BJTemitPrimeEmitPrimePtr) += (gpi + gepr + gm + go)
            //                                + (xcpi + xgm) * (s->real);
            // *(BJTemitPrimeEmitPrimePtr + 1) += (xcpi + xgm)
            //                                      * (s->imag);
            state.Matrix[BJTcolNode, BJTcolPrimeNode] += -gcpr;
            // *(BJTcolColPrimePtr) += (-gcpr);
            state.Matrix[BJTbaseNode, BJTbasePrimeNode] += -gx;
            // *(BJTbaseBasePrimePtr) += (-gx);
            state.Matrix[BJTemitNode, BJTemitPrimeNode] += -gepr;
            // *(BJTemitEmitPrimePtr) += (-gepr);
            state.Matrix[BJTcolNode, BJTcolPrimeNode] += -gcpr;
            // *(BJTcolPrimeColPtr) += (-gcpr);
            state.Matrix[BJTcolPrimeNode, BJTbasePrimeNode] += (-gmu + gm) + (-xcmu + xgm) * state.Omega;
            // *(BJTcolPrimeBasePrimePtr) += (-gmu + gm)
            //                                + (-xcmu + xgm) * (s->real);
            // *(BJTcolPrimeBasePrimePtr + 1) += (-xcmu + xgm)
            //                                     * (s->imag);
            state.Matrix[BJTcolPrimeNode, BJTemitPrimeNode] += (-gm - go) + (-xgm) * state.Omega;
            // *(BJTcolPrimeEmitPrimePtr) += (-gm - go)
            //                                  + (-xgm) * (s->real);
            // *(BJTcolPrimeEmitPrimePtr + 1) += (-xgm) *
            //                                      (s->imag);
            state.Matrix[BJTbasePrimeNode, BJTbaseNode] += (-gx);
            // *(BJTbasePrimeBasePtr) += (-gx);
            state.Matrix[BJTbasePrimeNode, BJTcolPrimeNode] += (-gmu) + (-xcmu - xcmcb) * state.Omega;
            // *(BJTbasePrimeColPrimePtr) += (-gmu)
            //                             + (-xcmu - xcmcb) * (s->real);
            // *(BJTbasePrimeColPrimePtr + 1) += (-xcmu - xcmcb)
            //                                     * (s->imag);
            state.Matrix[BJTbasePrimeNode, BJTemitPrimeNode] += (-gpi) + (-xcpi) * state.Omega;
            // *(BJTbasePrimeEmitPrimePtr) += (-gpi)
            //                                   + (-xcpi) * (s->real);
            // *(BJTbasePrimeEmitPrimePtr + 1) += (-xcpi)
            //                                      * (s->imag);
            state.Matrix[BJTemitPrimeNode, BJTemitNode] += (-gepr);
            // *(BJTemitPrimeEmitPtr) += (-gepr);
            state.Matrix[BJTemitPrimeNode, BJTcolPrimeNode] += (-go) + (xcmcb) * state.Omega;
            // *(BJTemitPrimeColPrimePtr) += (-go)
            //                                   + (xcmcb) * (s->real);
            // *(BJTemitPrimeColPrimePtr + 1) += (xcmcb)
            //                                    * (s->imag);
            state.Matrix[BJTemitPrimeNode, BJTbasePrimeNode] += (-gpi - gm) + (-xcpi - xgm - xcmcb) * state.Omega;
            // *(BJTemitPrimeBasePrimePtr) += (-gpi - gm)
            //                           + (-xcpi - xgm - xcmcb) * (s->real);
            // *(BJTemitPrimeBasePrimePtr + 1) += (-xcpi - xgm - xcmcb)
            //                  * (s->imag);
            state.Matrix[BJTsubstNode, BJTsubstNode] += (xccs) * state.Omega;
            // *(BJTsubstSubstPtr) += (xccs) * (s->real);
            // *(BJTsubstSubstPtr + 1) += (xccs) * (s->imag);
            state.Matrix[BJTcolPrimeNode, BJTsubstNode] += (-xccs) * state.Omega;
            // *(BJTcolPrimeSubstPtr) += (-xccs) * (s->real);
            // *(BJTcolPrimeSubstPtr + 1) += (-xccs) * (s->imag);
            state.Matrix[BJTsubstNode, BJTcolPrimeNode] += (-xccs) * state.Omega;
            // *(BJTsubstColPrimePtr) += (-xccs) * (s->real);
            // *(BJTsubstColPrimePtr + 1) += (-xccs) * (s->imag);
            state.Matrix[BJTbaseNode, BJTcolPrimeNode] += (-xcbx) * state.Omega;
            // *(BJTbaseColPrimePtr) += (-xcbx) * (s->real);
            // *(BJTbaseColPrimePtr + 1) += (-xcbx) * (s->imag);
            state.Matrix[BJTcolPrimeNode, BJTbaseNode] += (-xcbx) * state.Omega;
            // *(BJTcolPrimeBasePtr) += (-xcbx) * (s->real);
            // *(BJTcolPrimeBasePtr + 1) += (-xcbx) * (s->imag);
        }

        /// <summary>
        /// Load noise sources for the bipolar transistor
        /// </summary>
        /// <param name="ckt"></param>
        /// <param name="mode"></param>
        /// <param name="operation"></param>
        /// <param name="data"></param>
        /// <param name="OnDens"></param>
        /// <param name="p"></param>
        public override void NoiseLoad(Circuit ckt,
            Noise.OutputMode mode,
            Noise.Operation operation,
            NoiseData data,
            ref double OnDens,
            Model mod, object p = null)
        {
            BipolarModel model = mod as BipolarModel;

            CircuitStateComplex state = ckt.StateComplex;
            CircuitState cstate = ckt.State;
            NoiseJob curjob = ckt.CurrentJob as NoiseJob;
            double tmp = 0.0;
            double tempOnoise, tempInoise;
            double[] noizDens = new double[BJTnVar.Length];
            double[] lnNdens = new double[BJTnVar.Length];
            string[] BJTnNames = new string[] { "_rc", "_rb", "_re", "_ic", "_ib", "_1overf", "" };

            switch (operation)
            {

                case Noise.Operation.N_OPEN:

                    /* see if we have to to produce a summary report */
                    /* if so, name all the noise generators */

                    if (curjob.NStpsSm != 0)
                    {
                        switch (mode)
                        {

                            case Noise.OutputMode.N_DENS:
                                for (int i = 0; i < BJTnVar.Length; i++)
                                    data.namelist.Add(String.Format("onoise_{0}{1}", Name, BJTnNames[i]));
                                break;

                            case Noise.OutputMode.INT_NOIZ:
                                for (int i = 0; i < BJTnVar.Length; i++)
                                {
                                    data.namelist.Add(String.Format("onoise_total_{0}{1}", Name, BJTnNames[i]));
                                    data.namelist.Add(String.Format("inoise_total_{0}{1}", Name, BJTnNames[i]));
                                }
                                break;
                        }
                    }
                    break;

                case Noise.Operation.N_CALC:
                    switch (mode)
                    {

                        case Noise.OutputMode.N_DENS:
                            NevalSrc(ckt, ref noizDens[BJTRCNOIZ], ref lnNdens[BJTRCNOIZ],
                                     Noise.Sources.THERMNOISE, BJTcolPrimeNode, BJTcolNode,
                                     model.BJTcollectorConduct * BJTarea);

                            NevalSrc(ckt, ref noizDens[BJTRBNOIZ], ref lnNdens[BJTRBNOIZ],
                                     Noise.Sources.THERMNOISE, BJTbasePrimeNode, BJTbaseNode,
                                     cstate.States[0][BJTstate + BJTgx]);

                            NevalSrc(ckt, ref noizDens[BJT_RE_NOISE], ref lnNdens[BJT_RE_NOISE],
                                     Noise.Sources.THERMNOISE, BJTemitPrimeNode, BJTemitNode,
                                     model.BJTemitterConduct * BJTarea);

                            NevalSrc(ckt, ref noizDens[BJTICNOIZ], ref lnNdens[BJTICNOIZ],
                                     Noise.Sources.SHOTNOISE, BJTcolPrimeNode, BJTemitPrimeNode,
                                     cstate.States[0][BJTstate + BJTcc]);

                            NevalSrc(ckt, ref noizDens[BJTIBNOIZ], ref lnNdens[BJTIBNOIZ],
                                     Noise.Sources.SHOTNOISE, BJTbasePrimeNode, BJTemitPrimeNode,
                                     cstate.States[0][BJTstate + BJTcb]);

                            NevalSrc(ckt, ref noizDens[BJTFLNOIZ], ref tmp,
                                     Noise.Sources.N_GAIN, BJTbasePrimeNode, BJTemitPrimeNode,
                                     (double)0.0);
                            noizDens[BJTFLNOIZ] *= model.BJTfNcoef *
                                                   Math.Exp(model.BJTfNexp *
                                                       Math.Log(Math.Max(Math.Abs(cstate.States[0][BJTstate + BJTcb]), Noise.N_MINLOG))) /
                                                   data.freq;
                            lnNdens[BJTFLNOIZ] =
                                Math.Log(Math.Max(noizDens[BJTFLNOIZ], Noise.N_MINLOG));

                            noizDens[BJTTOTNOIZ] = noizDens[BJTRCNOIZ] +
                                                   noizDens[BJTRBNOIZ] +
                                                   noizDens[BJT_RE_NOISE] +
                                                   noizDens[BJTICNOIZ] +
                                                   noizDens[BJTIBNOIZ] +
                                                   noizDens[BJTFLNOIZ];
                            lnNdens[BJTTOTNOIZ] =
                                Math.Log(noizDens[BJTTOTNOIZ]);

                            OnDens += noizDens[BJTTOTNOIZ];

                            if (data.delFreq == 0.0)
                            {

                                /* if we haven't done any previous integration, we need to */
                                /* initialize our "history" variables                      */

                                for (int i = 0; i < BJTnVar.Length; i++)
                                {
                                    BJTnVar[i].LNLSTDENS = lnNdens[i];
                                }

                                /* clear out our integration variables if it's the first pass */

                                if (data.freq == curjob.NstartFreq)
                                {
                                    for (int i = 0; i < BJTnVar.Length; i++)
                                    {
                                        BJTnVar[i].OUTNOIZ = 0.0;
                                        BJTnVar[i].INNOIZ = 0.0;
                                    }
                                }
                            }
                            else
                            {   /* data.delFreq != 0.0 (we have to integrate) */

                                /* In order to get the best curve fit, we have to integrate each component separately */

                                for (int i = 0; i < BJTnVar.Length; i++)
                                {
                                    if (i != BJTTOTNOIZ)
                                    {
                                        tempOnoise = NoiseIntegration.Integrate(noizDens[i], lnNdens[i],
                                                                BJTnVar[i].LNLSTDENS, data);
                                        tempInoise = NoiseIntegration.Integrate(noizDens[i] * data.GainSqInv,
                                                                lnNdens[i] + data.lnGainInv,
                                                                BJTnVar[i].LNLSTDENS + data.lnGainInv,
                                                                data);
                                        BJTnVar[i].LNLSTDENS = lnNdens[i];
                                        data.outNoiz += tempOnoise;
                                        data.inNoise += tempInoise;
                                        if (curjob.NStpsSm != 0)
                                        {
                                            BJTnVar[i].OUTNOIZ += tempOnoise;
                                            BJTnVar[BJTTOTNOIZ].OUTNOIZ += tempOnoise;
                                            BJTnVar[i].INNOIZ += tempInoise;
                                            BJTnVar[BJTTOTNOIZ].INNOIZ += tempInoise;
                                        }
                                    }
                                }
                            }
                            if (data.prtSummary)
                            {
                                for (int i = 0; i < BJTnVar.Length; i++)
                                {     /* print a summary report */
                                    data.outpVector[data.outNumber++] = noizDens[i];
                                }
                            }
                            break;

                        case Noise.OutputMode.INT_NOIZ:        /* already calculated, just output */
                            if (curjob.NStpsSm != 0)
                            {
                                for (int i = 0; i < BJTnVar.Length; i++)
                                {
                                    data.outpVector[data.outNumber++] = BJTnVar[i].OUTNOIZ;
                                    data.outpVector[data.outNumber++] = BJTnVar[i].INNOIZ;
                                }
                            }    /* if */
                            break;
                    }    /* switch (mode) */
                    break;

                case Noise.Operation.N_CLOSE:
                    return;         /* do nothing, the main calling routine will close */
            }    /* switch (operation) */
        }
    }
}
